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The so-called Mixed Feedback Loop (MFL) is a small two-gene network where protein A regulates 
the transcription of protein B and the two proteins form a heterodimer. It has been found to be 
statistically over-represented in statistical analyses of gene and protein interaction databases and to 
lie at the core of several computer-generated genetic networks. Here, we propose and mathematically 
study a model of the MFL and show that, by itself, it can serve both as a bistable switch and as 
a clock (an oscillator) depending on kinetic parameters. The MFL phase diagram as well as a 
detailed description of the nonlinear oscillation regime are presented and some biological examples 
are discussed. The results emphasize the role of protein interactions in the function of genetic 
modules and the usefulness of modelling RNA dynamics explicitly. 

PACS numbers: 87.17.Aa, 87.16.Yc, 82.40.Bj 



I. INTRODUCTION 



Biological cells rely on complex networks of biochemi- 
cal interactions. Large scale statistical analyses have re- 
vealed that the transcriptional regulation networks and 
the networks of protein-protein interaction for different 
organisms are far from random and contain significantly 
recurring patterns pj or motifs. Mathematical modelling 
is useful to determine if these motifs can by themselves 
fulfill useful functions and it has, for instance, helped 
to show that the common 'Feed Forward Loop' mo- 
tif in transcriptional regulation can act as a persis- 
tence detector 0. More recently, a combined analysis of 
protein-protein interactions and transcriptional networks 
in yeast (Saccharomyces cerevisiae) has pointed out sev- 
eral motifs of mixed interactions H 0| . The simplest 
such motif composed of both transcriptional and protein- 
protein interactions is the two-protein Mixed Feedback 
Loop (MFL) depicted in Fig. ^ It is composed of a 
transcription factor A, produced from gene g a , and of 
another protein B , produced from gene gb- A regulates 
the transcription of gene gb and also directly interacts 
with protein B. The MFL has independently been ob- 
tained as the core motif of several networks produced in a 
computer evolutionary search aiming at designing small 
functional genetic modules performing specific functions 
Here, to better understand the possible functions of 
this basic module, we propose a model of the MFL based 
on the simplest biochemical interactions. This mathe- 
matical model is described in section [H] This is used to 
show that there exists wide ranges of kinetic parameters 
where the MFL behaves either as a bistable switch or 
as an oscillator. For the convenience of the reader, an 
overview of the different dynamical regimes is provided 



in section IIIII They are delimited in parameter space 
and their main characteristics are summarized. We then 
give detailed descriptions of the bistable regime in sec- 
tion [W] and of the nonlinear oscillations in section [V] A 
comparison is also made with a simple auto-inhibitory 
gene model with delay. In the concluding section, the 
important role played by protein dimerization and the 
usefulness of explicitly modelling mRNA dynamics are 
underlined and biological examples of the two proposed 
functions of the MFL are discussed. 



II. A MATHEMATICAL MODEL OF THE MFL. 
A. Model formulation 

As previously described, the MFL consists of two pro- 
teins A and B and their genes g a and gb, such that A 
regulates the transcription of gene gb and also directly 
interacts with B. Our aim is to analyze the dynamics of 
this small genetic module and see what can be achieved 
in the simplest setting. Therefore, different cellular com- 
partments and separate concentrations for the nucleus 
and cytoplasm are not considered and biochemical reac- 
tions are modelled by simple rate equations. 

The proposed MFL model is represented schematically 
in Fig.|3and consists of four equations that are described 
and explained below. The concentration of a chemical 
species X is denoted by square brackets and the cell vol- 
ume is taken as volume unit. So, [X] represents the ef- 
fective number of X molecules present in the cell. 

The first two equations model the transcriptional reg- 
ulation of gene gb by protein A: 
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6[g b :A]-a[g b }[A] (1) 
Pf[9b] + Pb[gb ■ A] - S T [i h ] (2) 
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where it is assumed that gene gb exists under two forms, 
with A bound to its promoter with probability [g b : A] 
and without A with probability [gb]- Since [gb] + [gb ■ 
A] = 1, the single Eq. is sufficient to describe the 
transition between the two forms. Specifically, A proteins 
bind to the gb promoter at a rate a and when bound they 
are released at a rate 9. The regulation of transcription 
of gene gb by protein A is described by Eq. J2J where 
rb stands for gb transcripts. When A is bound to gb pro- 
moter, transcription is initiated at a rate pb and otherwise 
it is initiated at a rate pf. Thus, pb > Pf corresponds 
to transcriptional activation by A and pb < Pf to tran- 
scriptional repression. A first order degradation for gb 
mRNA at a constant rate 5 Y has also been assumed. As 
given, the description is strictly valid for a single copy 
gene. However, it also applies for a gene with go copy 
(e.g go — 2) provided that pf and pb are understood 
to be go-fold greater than the corresponding elementary 
rates. 

The production of A and B proteins and their com- 
plexation are described by the following two equations 
that complete our description of the MFL module, 



dJA] 
dt 

d[B] 
dt 



PA-7P] [A]-* A [A] 
+6[g b : A]-a[g b ][A] 

/3[r b ]- 7 [B] [A]-J B [B] 



(3) 
(4) 



where [A] and [B] respectively denote the concentration 
of proteins A and B. Since, regulation of gene g a is not 
considered, separate descriptions of its transcription and 
translation are not needed and it is simply assumed in 
Eq. © that protein A is produced at a given basal rate 
PA- The second crucial interaction of the MFL, the direct 
interaction between protein A and B is taken into account 
by assuming that A and B associate at a rate 7. It is as- 
sumed that B does not interact with a protein A that 
is bound to gb promoter. In addition Eq. assumes a 
first-order degradation of protein A at a constant rate 5 a, 
and its last two terms come from the (small) contribu- 
tion to the concentration of A in solution of the binding 
(unbinding) of A to (from) g b promoter. Eq. assumes 
that B proteins are produced from the transcripts of gene 
gb at a rate j3 and are degraded at a rate 5b in addition 
to their association with A proteins. 

As described by Eq. l.'il l[l . the complexation between 
A and B proteins proceeds at a rate 7. For simplicity, 
we suppose that the AB complex does not interact with 
genes g a and gb, we neglect its dissociation [f| and simply 
assume that it is fully degraded at a rate 5ab after its 
formation, 



d[AB] 
dt 



= 7 [B] [A]-6 AB [AB] 



(5) 



Since the complex AB does not feed back on the dynamics 
of the other species, its concentration does not need to 
be monitored and Eq. JSJ is not explicitly considered in 
the following. 



B. Values of kinetic parameters and a small 
dimensionless parameter 

Even in this simple model, ten kinetic constants should 
be specified. It is useful to consider the possible range 
of their values both to assess the biological relevance of 
the different dynamical regimes and to orient the model 
analysis. 

Half- lives of mRNA range from a few minutes to several 
hours and are peaked around 20 minutes in yeast 0. 
S r = 0.03m?n _1 can therefore be taken as a typical value. 

For the transcription factor-gene promoter interaction, 
typical values appear to be a critical concentration 9/a = 
[A]q in the nanomolar range, a bound state lifetime of 
several minutes and activated transcription rates of a few 
mRNAs per minute. We therefore assume a = 0/40, that 
9 is of the same order as 6 r and pf ~ min -1 , pb ~ min^ 1 . 

Proteins half-lives vary from a few minutes to several 
days The hour appears as a typical value. We choose 
Sa = 6b = O.Olmin -1 and more generally consider that 
the A and B protein half-lives are of the same order or 
longer than that of gb mRNA. For protein production, 
j3 = 3 protein molecules per mRNA molecule per minute 
appears a plausible value for an eucaryotic cell 9] . This 
gives pa — 100mm -1 when combined with the previous 
values for mRNA production. 

We assume that protein interaction is essentially lim- 
ited by diffusion. A diffusion constant D ~ 2.5fj,m 2 s~ 1 . 
[iril | gives a time s 2 / D of about a minute for diffusion 
across a cell of a size s — 10 p,m. We therefore choose 
7 ~ I min . 

It is convenient to adimension Eqs. 11 l.'il) to decrease as 
far as possible the number of independent parameters. 

We first normalize the g b mRNA concentration by the 
concentration that gives a production of B protein equal 
to that of A. We thus define the dimensionless concen- 
tration r — /3[rb]/pA- We normalize the protein concen- 
trations by the equilibrium concentration resulting from 
production at a rate pa and dimerization at a rate 7 and 
define A = ^/7/pA [A], B = \J^j p A [B]. We also take as 
time unit the inverse of gb mRNA degradation rate l/5 r . 
With theses substitutions, Eqs. Ill II) read, 
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where we have defined the following rescaled parame- 
ters 6 = 5 r /^/p A j , pi = (3p b /(p A 5 r ), po = Pp f /(p A 6 r ) 

9 = 9/6 r , /i = y/i/pA, d a = 5 A /6 r and d b = 6 B /6 r . 
We have also introduced th e dime nsionless critical bind- 

\JllPA 0/a 



ing concentration Aq 
depends on seven parameters. 



The model still 
In order to simplify its 
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analysis, it is useful to note that S is a small parameter 
(approximately equal to 3 x 10~ 3 with the previous esti- 
mations) . The influence of three key parameters of order 
one in Eq. 1GI9I) is particularly examined in the following. 
These are po and pi which measure the strengths of the 
two possible states of B production (with or without A 
bound to gene gb) as compared to that of A, and 9 which 
compares the rates of A unbinding from DNA to that 
of mRNA degradation (a is supposed to vary with 9 to 
maintain a fixed critical binding concentration Aq). 



III. OVERVIEW OF THE DYNAMICS IN 
DIFFERENT PARAMETER REGIMES. 

We provide here an overview of the different dynamical 
regimes of the MFL which are depicted in Fig.[3]and sum- 
marize their characteristics. The behavior of the MFL 
depends on whether protein A is a transcriptional acti- 
vator or a transcriptional repressor and on the strengths 
of the two transcription rates of gene gb (i.e. with A 
bound or not to its promoter). More precisely, the key 
parameters are the strengths of B protein production, 
j3pf/5 r and /3pb/6 r , as compared to the production rate 
Pa of protein A. It is therefore convenient to consider 
the previously introduced ratios po — (3pf /(S r pA) and 
Pi = f3p b /(S r p A ). 

We have observed that depending on the values of po 
and p\ the MFL can be monostable, exhibit bistability 
or display oscillations. We qualitatively describe these 
three different cases in the following. Their biological 
relevance is further discussed in section IVTl 



A. Monostable steady states 

The simplest case occurs when the production rate of 
A is higher or lower than the production rate of B, irre- 
spective of the state of gb promoter. That is when both 
/3pf/S r and [3pb/8 r are either higher or lower than pa- In 
this case, the MFL has a single stable state to which it 
relaxes starting from any initial conditions. 

When both B production rates are higher than the pro- 
duction rate of A (i.e. po > 1 and p\ > 1), A proteins are 
quickly complexed by B and are unable to interact with 
gb promoter. The concentration [A] of uncomplexed A 
proteins is therefore low and results from a simple bal- 
ance between production and complexation. The high 
concentration of uncomplexed B proteins is the effective 
result from transcription at the free gb promoter rate, 
complexation and degradation, 




An equally simple but opposite result holds when both 
B production rates are smaller than the production rate 
of A (po < 1 and p\ < 1). Then, the concentration of 



uncomplexed A is high , gb promoter is occupied by A 
and a low B concentration results from a balance between 
complexation and degradation 




The dynamics of the MFL is richer when the produc- 
tion rate of A is intermediate between the two possible 
production rates of B, [3pf/8 r and f3pb/S r i.e. when ei- 
ther po > 1 > Pi or pi > 1 > po. We consider these two 
cases in turn. 



B. Transcriptional repression and bistability. 

We first discuss the case when A is a transcriptional 
repressor, pQ > 1 > pi. Then, two stable steady states 
can coexist. Let us suppose first that no A is bound to gb 
promoter. Then the production rate of B is larger than 
the production rate of A, and all produced A proteins are 
quickly complexed. This stably prevents the binding of A 
proteins to gb promoter and maintain a steady state with 
low A and high B concentrations approximately equal to 

M'" 5- [A] ^Ww' (12) 

The second opposite possibility is that A is sufficiently 
abundant to repress the transcription of gene gb- Then, 
since the production rate of A has been supposed to be 
higher than the production rate of B in the repressed 
state, B proteins are quickly complexed but uncomplexed 
A proteins are present to maintain the repression of gene 
gb transcription. This gives rise to a second stable state 
with high A and low B concentrations approximately 
equal to 

The bistability domain is only exactly given by the 
simple inequalities pq > 1 and pi < 1 when the ratio S 
of gb mRNA degradation rate over the effective protein 
dynamics rate is vanishingly small. As shown in section 
l^and on Fig.|3 for small value of S, it is more accurately 
given by 

Po > 1 + 2^(1- Pl )S B /A Q , Pl < 1, 

p x < 1 - 2y/(pa-l)Ao6A, Po > 1- (14) 

In this intermediate range of production of A proteins, 
the network reaches one or the other fixed points, de- 
pending on its initial condition. The existence of this 
bistability domain can serve to convert a graded increase 
in A production in a much more abrupt switch-like re- 
sponse in A (and B) concentration as shown in Fig. 01 
The usefulness of this general feature of multistability 
has been recently discussed in different contexts [ill Il2| . 
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C. Transcriptional activation and oscillations. 

When A is a transcriptional activator, the complexa- 
tion of B with A acts as a negative feedback and can 
serve to diminish the variation in B protein concentra- 
tion when A varies. This leads also to oscillations when 
A production lies in the intermediate range p\ > 1 > po. 
This oscillatory behavior mainly depends on the ratios 
of protein production po, p\ and exists for a large range 
of DNA-protein interaction kinetics . However, a faster 
kinetics leads to a smaller oscillatory domain as shown 
in Fig. Oscillations cannot be sustained when 8 be- 
comes large and comparable to 1/5 (with Aq fixed), or 
equivalently when 9 ~ S r /S. It can also be noted that 
for a sufficiently large activation of transcription by A, 
there exists a domain of coexistence between oscillations 
and steady behavior as shown in Fig. he. depending 
on initial conditions concentrations will oscillate in time 
or reach steady levels. 

For small 5, oscillations are nonlinear for most param- 
eters. As can be seen in Fig. El an oscillation cycle com- 
prises two phases in succession, a first phase of duration 
Ti when the concentration of protein A is high followed 
by a second phase of duration Ti where the concentration 
of B is high. A full oscillation cycle of period T = T\ +T2 
proceeds as follows. Let us start with low concentrations 
of A and B proteins at the beginning of the first phase. 
When no A is bound to gb promoter, B production rate 
is lower than A production rate and complexation can- 
not prevent the increase of A concentration. When the 
concentration of A has reached a critical level, A start to 
bind to gb promoter and activate transcription. This re- 
sults in a higher production of B than A and the diminu- 
tion of free A by complexation. Since A concentration 
is high, the produced B's are quickly complexed and the 
concentration of uncomplexed B's remain low. Eventu- 
ally, A concentration drops below the binding level and 
no longer activates gb transcription. This marks the tran- 
sition between the two parts of the oscillation cycle. B 
continues for a while to be produced from the transcripts 
and since few A's are present this now leads to a rise of 
the concentration of free B's. Finally, the concentrations 
of B transcripts and B proteins drop and phase II of the 
oscillation cycle terminates. A new cycle begins with low 
concentrations of A and B proteins. 

One remarkable feature of the nonlinear oscillations 
coming from the smallness of the parameter 5, is that the 
concentrations of A and B proteins in their respective 
high phase depend weakly on the complex association 
rate 7, as long as it is not too small for the oscillation to 
exist, and that the period of the oscillations is strikingly 
insensitive to the exact value of 7. For the parameters 
of Fig. [f)| the oscillation period only changes by about 
1% when 7 is reduced from lmin^ 1 to 2 x 10~ 2 min . 
This remains true even if the formation of the complex 
AB is not irreversible as supposed in the present model. 
We have checked that the case when the complex AB 
forms with an association rate 7 Q , dissociates with a rate 



-fd and is degraded with a rate Sab is well described 
by the present model when one takes the effective rate 
7 = 7<j#ab/(7c2 + Sab) for the irreversible formation of 
the complex (as obtained from a quasi-equilibrium as- 
sumption). 

In the two following sections, we provide a more de- 
tailed analysis of the MFL different dynamical regimes 
and derivations of the results here summarized. 



IV. THE SWITCH REGIME. 

We first determine the possible steady states for the 
MFL. The free gene, mRNA and B protein concentra- 
tions are given in a steady state as a function of A con- 
centration by, 



9 = 



B 



A 



A + A Q 
piA + pqA 
A + A 
PiA + p A 
(A + A Q )(A + Sd b ) 



(15) 
(16) 
(17) 



The concentration of A itself satisfies the following equa- 
tion 



1 = Sd a A 



( Pl A + p A )A 
(A + A Q )(A + Sd b ) 



(18) 



For 5 — 0, the right-hand side (r. h. s. ) goes monotoni- 
cally from po at small A to pi at large A and a solution, 
A2 = Aq(pq — 1)/(1 — pi) exists only when po < 1 < pi 
or pi < 1 < po- For small 5, two other steady states are 
possible. A steady state with a small concentration of A, 
Ai ~ 5db/(po — 1), exists when p > 1. Inversely a steady 
state with a large concentration of A, A3 ~ (1 — p\)/{Sd a ) 
exists when p\ < 1. Therefore, Eq. I|18fl has multiple (i.e. 
three) fixed points only when A is a transcriptional re- 
pressor in the region p\ < 1 < pq. This is the parameter 
regime which we examine in the rest of this section. 



A. The high A state 

We show that the state with a high concentration of A 
proteins is stable. The form of Eqs. © and © suggests 
to simplify the analysis by using the fact that protein 
quickly reach a quasi-equilibrium state. However, both 
proteins cannot be in quasi-equilibrium at the same time. 
For instance, when proteins A 3> B, only B is in quasi- 
equilibrium and A scale as 1/5. When A concentration is 
high, in the limit of small 5, we set a = 5 A. Substitution 
in Eq. (Ill II) shows that g and B reach on a fast time-scale 
their quasi-equilibrium concentration 



9 
B 



5A /a 
5r/a 



(19) 
(20) 
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Therefore the dynamics of the MFL reduces to the fol- 
lowing two equations in the limit S — ► : 



di- 
di 
da 

dt 



pi-r 

1 — r — d„a 



(21) 
(22) 



Eqs. i|2H22l) clearly show that the high A fixed point 
is stable and as found above satisfies r ~ pi, B ~ g ~ 0, 
a = (1 — p\)/d a at the lowest order. This steady state 
is possible only if pi < 1, that is when B production 
rate is not high enough to titrate A and to prevent the 
repression by A. 



B. The high B state. 

The high B state can be analyzed in a very similar way 
in the limit of small 5. When B concentration is high, 
we define b = 5B. In that case, A quickly reaches its 
quasi-equilibrium concentration: 



A ~5/b 



(23) 



and the MFL dynamics reduces to the following three 
equations: 



dg 
dt 
dr 

dt 
db 

dt 



0(1-9) 

Po5 + Pl(l - 9) -r 
r — 1 — di,b 



(24) 
(25) 
(26) 



The concentrations tend toward those of the high B 
fixed point g ~ 1, r ~ po, b ~ (pa — l)/d& A ~ 0. This 
steady state exists only if po > 1, when the production 
of B proteins is high enough to titrate A proteins and to 
prevent transcriptional repression by A. 



Bifurcation between the monostable and 
bistable regimes. 



For A ~ \f5 Eq. (|18|l becomes at dominant order after 
expansion, 



. A db 

p °- 1 = (1 - pi) T a +5 A 



(27) 



For a given repressio n strength pi < 1, the r. h. s. has a 
minimum for A = \J SdbAoJJl — pi) (that is of order ^/S 
as assumed) such that 



Po 



ISd b (l-pi) 



Ao 



0(5) 



(28) 



Above this value of po, Eq. I|27(l has two solutions and 
below it has none. Eq. i|28|l therefore marks the saddle- 
node bifurcation line that separates the bistable domain 
from the monostable domain with a high A concentration 
fixed point. For S — ► 0, the zeroth order boundary po = I 
is of course retrieved. 

The approximate expression (|28|l agrees well with an 
exact numerical determination of the bifurcation line as 
shown in Figure 

The transition between the bistable regime and the low 
A concentration monostable regime (top-right corner of 
Fig. |3J) can be similarly analyzed. The high (A3) and 
middle (A2) A concentration fixed points can coalesce 
when for both the concentration A is of order 1 / y/6. Ex- 
pansion of Eq. (|18|) for A ~ 1/ y/S gives 



1 - Pi = (po 



l)^- + 6d A A 



(29) 



For a given p > 1, the r. h. s. of Eq. I|29|) is minimum 
for A = y/ Aq(pq — l)/(8d a ) and at this minimum pi is 
equal to 



Pi 



= l-2y/d a Ao6(pa-l) 



(30) 



For pi smaller than this value, Eq. I|29|l has two roots 
corresponding to A2 and A3. The two roots merge and 
disappear on the saddle-node bifurcation line (|3L)fl that 
marks the boundary between the bistable domain and 
the low A concentration monostable regime. Comparison 
between the approximate expression l|30(l and an exact 
numerical determination of the bifurcation line is shown 
in Figure |3| 



The previous analysis has shown that bistability exists 
in the domain po > 1 > Pi when i « 1. We analyze 
more precisely the nature and position of the bifurcations 
between the monostable and bistable regimes. 

We consider first the transition between the bistable 
region and the monostable region with a high concen- 
tration of A (bottom left corner in Fig. |2J. It follows 
from Eq. JTHJ) that the low (A x ) and middle (A 2 ) A 
concentration fixed points coalesce and disappear when 
po — 1 = 0(V~5) and both A± and A2 are of order \~8. 



V. THE OSCILLATOR REGIME 

We consider now the case when A is a transcriptional 
activator, that is when po < p\. In this case, the MFL 
has a single steady state since the right-hand side (r.h.s.) 
of Eq. I|18l) is a monotonically increasing function. How- 
ever, as we show below, this steady state is unstable in a 
large domain of parameters and the MFL behaves as an 
oscillator. 
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A. The linear oscillatory instability 

We begin by analyzing the linear stability of the fixed 
point (g*, r*, B* , A*) where A* is the single solution of 
Eq. H18fl and g*,r* and B* are given as functions of A* 



by Eqs. (|f 5117(1 . After linearization of the MFL dynamics 
[Eqs. I(6I9( around this fixed point, the complex growth 
rates a of perturbation growing (or decreasing) exponen- 
tially in time like exp(cri) are found to be the roots of the 
following characteristic polynomial, 



A*. 
A~v 



B* 

T 



A* 

T 



\ A*B*~ 


g*9 


) S 2 


"AT 



r 



A* ( A* \ 

£2 (A) -Pi) -M<7 {a + d b + — J 

(311 



Again, the fact that 5 is small simplifies the analysis. 
We consider first the case when the rescaled concentra- 
tions g* ,r* , B* , A* are of order one. This corresponds to 
the fixed point A2 of the previous part in the parameter 
regime p < 1 < pi, with A* ~ A (l - p )/(pi - 1), 
B* ~ 1/A* and g* ~ (pi — l)/(pi ~ Po)- Let us assume 
that the roots cr diverge as S tends to zero. Then, Eq. I(31|) 
reduces at dominant order to 



a 3 (S<T + A*+B*) = ^^(p a - Pl ) (32) 
A 

Therefore, three roots Uk, k — 0,1, 2, are of order 5 -1 / 3 
and proportional to the three cubic roots j k of unity 



07c 



.4* 



A* + B 



9' 0, 

T 8^ 



Po) 



1/3 



k = 0,1,2. 



The fourth root <t 3 is of order 1/5 and corresponds to a 
stable mode of real part —(A* + B*)/5. The two roots 
1 and 02 are complex conjugate and have a positive real 
part. Thus, in this parameter domain, the MFL fixed 
point is oscillatory instable. The dynamics tends toward 
an attractive limit cycle that we will analyze in the next 
subsection. When A* or B* grows (i.e when p\ — ► 1 or 
po — ► 1) the fixed point instability disappears via a Hopf 
bifurcation. We analyze these two boundaries in turn. 



positive real parts. On the contrary, when pi becomes 
sufficiently close to one, the r.h.s of Eq. ((34(1 becomes 
negligible and Eq. (|3"P has obviously three real negative 
roots. Therefore, when p\ — > 1, one traverses the bound- 
ary of the linearly unstable region. Two roots of Eq. 1(34|) 
traverses the imaginary axis on this boundary in a Hopf 
bifurcation. Assuming (and checking afterwards) that, 
on this boundary, these roots are small compared with 
l/(Pi — 1)> their expressions is found perturbatively by 
expanding the first factor in Eq. I|34|). 



i + d a 



(pi - l) s 



_ ± i / (pi-i) 2 



2 2(pi - po) 2 A 9S 2\A {p 1 -po)5 

(35) 

This provides the location of the stability boundary, 



Pi = 1 + S9A (d a + l)(pi - po)' 



1/3 



(36) 



In the limit <5 — s- the zeroth order boundary p\ = 1 is 
of course retrieved. It can be checked that the obtained 
result pi — 1 ~ t) 1 / 3 justifies a posteriori the use of the A2 
5- independent expression for the fixed point. This fixed 
point linear stability boundary 1)36(1 can also be directly 
obtained from the Routh-Hurwitz condition 01 on the 
third degree polynomial ((34(1 . Comparison of the small 
8 Eq. ((36(1 with the exact linear stability boundary is 
provided in Fig. |3J 



A* high 



2. 



A* small 



When pi — * 1, A* becomes large and the real parts of 
the roots become of order one. Thus, Eq. l(3*Tl) needs to 
be approximated in a different way. We neglect 6(a + df,) 
and B* as compared to A* in Eq. ((31|) and obtain 



(a + 



Pi ~ PQW . ,W r T \ ^Pl - 1) , q/ A 

)(<7 + l)(o- + d a ) = — (34) 



Pi 



1 



5A 



where we have replaced A* and g* by their expressions at 
the A 2 fixed point. When the r.h.s is the dominant con- 
stant in Eq. ((34(1 , its three roots are proportional as above 
to the three roots of (minus) unity and two of them have 



We now consider the upper boundary (po ~ 1) of the 
linearly unstable region. When po tends toward one, B* 
grows l/(po — 1) and becomes large compared to 8(a+d a ) 
and A*. Eq. I(31|) then simplifies and reduces at leading 
order to: 



with 



{a + 6)((j + l){a + d b ) = -E 



A B*S 



(37) 



(38) 



7 



Again when E is large the roots are proportional to the 
three roots of -1 and two have positive real parts. On 
the contrary, the three roots are real and negative when 
E vanishes. There is a critical value E c such that for 
E < E c the fixed point is linearly stable. This occurs 
when the real part of the two complex conjugate roots 
becomes negative at the value E c given by the Routh- 
Hurwitz criterion [T^j 



E c = (d b + l + 0)(§ + d b + d b 0) - 



(39) 



The allied value of A at the fixed point concentration is 
obtained from Eq. 



.4* 



' E C A S 
0(pi - 1) 



(40) 



Using the fixed point Eq. (|18fl . this in turn corresponds 
to the line in parameter space 



Po = 1 - 



'E c S(jn-l) 



And 



1 -i 



E c 



(41) 



Eq. Ij41|l , of course, reduces to the zeroth order boundary 
Po = 1 when 5 — > 0. Comparison of the small-5 asymp- 
totic expression (|41|l with the exact linear stability line 
is provided in Fig. |31 



B. A description of the nonlinear oscillations 

The MFL oscillations quickly become strongly nonlin- 
ear away from threshold (or, of course, when the bifur- 
cation is subcritical). As there is no auto regulatory di- 
rect or indirect positive feedback loops in the MFL, two- 
variable reductions have a negative divergence. It follows 
from the well-known Bendixson's criterion |l3j |. that they 
cannot be used to describe the oscillatory regime. It is 
for instance not possible to properly describe the oscilla- 
tions by only focusing on the protein dynamics. A spe- 
cific analysis is therefore required and developed in this 
section. 

A full period of the nonlinear oscillations can be split 
into two parts: a first phase with A » 1, and a second 
phase with B » 1 (see Fig. [7J|. 

In the following, these two main phases of the limit 
cycle are described. We use simple continuity arguments 
to match the two phases and obtain a description of the 
whole limit cycle. The oscillation period is computed (at 
the lowest-order in $). A full justification of this simple 
matching procedure and a detailed study of the transition 
regimes between the two main phases are provided in 
Appendix fusing matched asymptotics techniques. 



1. Phase I: High A, Low B phase. 

This phase is defined as the fraction of the limit-cycle 
where the concentration of protein A is larger than [A]q 



i. e. when A » 1 >> B. Eq. © leads us to assume 
that A is of the order 1/5. So we define a = 5 A. In 
the limit of small 5, as A scales as 1/5, both the binding 
to g promoter and the dynamics of B resulting from its 
complexation with A are fast compared to that of mRNA 
and A. Consequently, g and B are in quasi-equilibrium 
and obey at lowest order in 5, 



fJ 
B 



5A Q /a 
5r/a 



(42) 
(43) 



The dynamics of the MFL therefore reduces to the fol- 
lowing two equations at lowest order in 5, 



dr 

dt 
da 



Pi-r 

1 — r — d a a 



(44) 
(45) 



The beginning of phase I coincides with the end of phase 
II where, as explained below, A concentration is small. 
Continuity therefore requires that a vanishes at the start 
of phase I (a detailed study of the transition region be- 
tween the two phases is provided in Appendix [SJ . De- 
noting by ri, the value of r at the start of phase I, an 
easy integration of the linear Eq. 1|44I45(1 gives, 



ri(t) = Pi + (n - pi)e * 
1 

ai(t) = - 



n - pi 



(46) 
3-*] (47) 



Subscripts / have been added to r and a in Eq. i|4l)|) and 
(|47|l to emphasize that the corresponding expressions are 
valid during phase I only. Phase I ends with the fall of A 
concentration, after a time t\ such that ai(t\) = 0. The 
mRNA concentration is then equal to ri(ti) = r 2 . 

One can note that the rise and fall of the concentration 
of protein A imposes restrictions on the parameters. The 
rise at the beginning of phase I is possible only if A pro- 
duction dominates over its complexation with B, that is if 
r\ < 1 (Eq. (|45|l . The following fall requires the reverse 
which can only happen if B production becomes suffi- 
ciently important, namely r > 1. This requires r 2 > 1 
and a fortiori p\ > 1 (Eq. lEH ). 



2. Phase II : High B phase. 

This second phase is defined as the part of the limit- 
cycle where the concentration of protein A falls below 
[A] i.e. B » 1 » A. The form of Eq. (jSJ) leads us to 
assume that B scales as 1/5 and we define b = SB. 

In the limit of small 5, the dynamics of A is fast com- 
pared to that of the other species and A is in quasi- 
equilibrium. At lowest order its concentration reads, 



A ~5/b 



(48) 
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Thus, the dynamics of the MFL reduces, at lowest order 
in S, to the following three equations, 



dt 

dr 
dt 
db 

d7 



0(1 - 9) 

Po9 + Pi(l-g) ~r 
r — 1 — djjb 



(49) 
(50) 
(51) 



Continuity with the previous phase I leads us to require 
that at the beginning of phase II b — 0, r — r 2 and g = 
(see Appendix 1X1 for a detailed justification). With these 
boundary conditions, the linear Eqs. 0§l-|2] can readily 
be integrated to obtain: 



9i i (t) = 
ru(t) = 

bn(t) = 



1 - e 
Po + 
Po - 1 



-01 



7'2 - PO + 



Pi - PO 



di, 

r2 - Po + 



1 _ e -<V] _ 
Pi - Po 



1 



pi - po e~ dht - er et 



e-d b 



i 



e -d b t _ e -t 



1 - 



(52) 



Since b vanishes at the beginning of phase II, b should 
start by rising. This imposes that B production dom- 
inates over complexation and requires that the concen- 
tration r at the beginning of phase II is greater than 1, 
i. e. r 2 > 1 (see Eq. (15 1|) ^) . Thus, r always remains larger 
than po and would decrease toward po for a long enough 
phase II (Eq. lj5U|l and j53)). At the end of phase II, b 
should decrease to to continuously match the beginning 
of phase I. This requires that complexation then domi- 
nates over production of B that is r < 1 (Eq. (ED), and 
it is only possible when po < 1. With these conditions 
met, phase II lasts a time i 2 and ends when 6// (£2) = 0. 



3. Matching of the two phases and period determination 

In order to complete the description of the limit cycle, 
it remains to determine the four unknowns ri,r 2 ,ti and 
i 2 from the four conditions coming from the continuity 
of proteins and mRNA concentrations, 



ri(h) = r 2 , oj(ti) = 0, 
ru(t2) = n, bn(t 2 ) = 0. 



(53) 
(54) 



One possibility to solve these equations is to use Eq. 153|) 
to express r± and r 2 as a function of t\. It is then not 
difficult to see that the third equation l|54|) implicitly de- 
termines ti as a function of t\. Once ri, r 2 and i 2 are 
obtained as functions of t\, the last equation 6(£ 2 ) = 
can be solved for t\ by a one-dimensional root finding 
algorithm. In this way, we have determined r% , r 2 , t\ , t% 
for different sets of kinetic constants and obtained the 
rescaled period of oscillation 



T r — t\ + t 2 



(55) 



the dimensionfull period being T = T r /8 r . Results are 
shown in Fig. |SJ The analytical period compares well to 
results obtained by direct numerical integration of Eq. JJJ 
01 for different values of 5 with, of course, a closer agree- 
ment for smaller S. 

Eqs. 1)531541) are difficult to solve analytically but an- 
alytic expressions can be obtained for various limiting 
cases. For instance, if the degradation of protein B is 
negligible (db = 0), Eq. (|52*|l simplifies to 



bu(t) - (po-l)t- 



Pi - Po 



r-2 - Po + 



Pi 



-1) 
Po 



1 



(1 



-9t\ 



(56) 



If we further consider the limit of large p\ , the condition 
bufa) = 0, gives the following estimate for the duration 
of phase II, up to exponentially small terms, 



t-2 



T2 ~ PO Pi ~ PO 
1 - PO 0(1 - po) 



(57) 



as well as the transcript concentration at the end of phase 
II 



n ^ po 



(58) 



Similarly, the condition ai{t\) = together with Eq. I|47() 
show that t\ is small for large pi, 



t x ~ 2(1 - n)/ Pl ~ 2(1 - p )/pi 



(59) 



where expression (|58|) has been used in the second equal- 
ity. Given the duration l|59() of phase I, the concentra- 
tion r 2 of transcript at its end is directly obtained from 
Eq. g|| and (|EEJ 



''2 



PO 



(60) 



This finally provides the estimate of the period for large 
pi (with db — 0), 



to 



Pi - Po 
6(1 - Po) 



(61) 



A comparison between Eq. (|61|l and numerically deter- 
mined oscillation periods is provided in Fig. |SJ 



Comparison with a simple delayed negative 
feedback. 



Popular models for genetic oscillators consist in a pro- 
tein repressin g it s own production with a phcnomcnolog- 
ical delay [Hill 03 • 

Using a Hill- function to model this 
repression, the simplest model of this kind reads, 



dB(t) 
dt 



1 + [B(t - t)/BoY 



cB 



(62) 
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where p is the protein production rate, c the protein 
degradation rate and r the phenomenological delay for 
repression. The phase diagram of this simple oscilla- 
tor can be computed [l^. In the limit of long delays 
(ct 1), the oscillations are nonlinear. In an expansion 
in the small parameter e = 1/ cr, their period T is 

T/t = 2 + — + ■ ■ • (63) 

CT 

where k is a constant which is approximately equal to 2 
for n = 2 and p/c of order one [19|. 

It is interesting to note some analogies between this 
simple model with delay and the previously studied MFL. 
When t >> 1/c, the period of the delayed model scales 
as r while the period of the MFL scales as S^ 1 . The 
mRNA half-life in the MFL thus plays the role of the 
phenomenological delay in Eq. (|62|) . This is in line with 
the mRNA being the major pacemaker of the MFL os- 
cillator . The rescaled parameter e in the simple model 
with delay plays a role similar to the rescaled parameter 
S in the MFL in the sense that, first, for high values of 
e, oscillations disappear jl8| . and for small values of this 
parameter, the period of the oscillations is independent 
of 1/ct at dominant order. In the delayed model, e rep- 
resents the ratio between the typical life-time of protein 
over the delay, while in the MFL, 5 represents the ratio 
between the typical time scale of protein production and 
sequestration over the typical life-time of RNA. These 
analogies show that for the MFL JpAl plays the role of 
the protein degradation constant c in the delayed model. 
Thus protein sequestration by complexation is the MFL 
analog of protein degradation in the simple model with 
delay. The fastness of both these processes relative to 
a slow mechanism (delay r, RNA dynamics) ensures the 
flipping of both oscillators between two states (two con- 
centrations of B for the delayed model, A high and B 
high for the MFL). 

VI. DISCUSSION 

We have proposed a simple model of the Mixed 
Feedback Loop, an over-represented motif in different 
genomes. We have shown that by itself this motif can 
serve as a bistable switch or generate oscillations. 

A. Importance of dimerization and of RNA 
dynamics 

The MFL dynamics crucially depends on the post- 
transcriptional interaction between A and B proteins. 
The dynamics of dimerization is fast and allows a high 
concentration of only one of the two proteins, thus ef- 
fectively creating a dynamical switch between the two 
species. If, for instance, A proteins are present at a higher 
concentration than B proteins, all B's are quickly titrated 
and there only remains free A and AB dimers. 



A bistable system is obtained by coupling this post- 
transcriptional mechanism to a transcriptional repres- 
sion. The created bistable switch is quite different from 
the classical 'toggle switch' [2(J which is based on recipro- 
cal transcriptional repressions between two genes. In the 
MFL model, the single transcriptional repression of gene 
gb is sufficient for a working switch. Moreover, dimer- 
ization bypasses the need for cooperativity in the toggle 
switch 21], and may render the present module simpler 
to implement in a biological system. 

When the dimerization between A and B is coupled 
to a transcriptional activation, the system behaves as an 
oscillator if the production of A proteins is intermediate 
between the two possible for B. Then, the level of gb tran- 
scripts controls A protein dimerization and, of course, B 
production rates. At the beginning of the derepression 
phase (phase I), the low level of transcripts leads to a 
small production of B proteins and to a rise of A con- 
centration. This first phase ends when the subsequent 
high production and accumulation of gb transcripts give 
rise to a production of B proteins sufficient for the titra- 
tion of all free A's. During the ensuing repression phase 
(phase II) , gb transcription rate is low and the transcripts 
degradation produces a continuous decrease of their con- 
centration. The parallel decrease of B production rate 
ultimately leads to the end of this second phase when B 
production rate is no longer sufficient for the titration of 
all produced free A's. 

In summary, the core mechanism of the clock is built by 
coupling the rapid switch at the protein level to the slow 
mRNA dynamics. The protein switch in turn controls 
mRNA dynamics via transcriptional activation. 

It is interesting to note that the MFL model in the 
simple form here analyzed provides a genetic network 
that oscillates without the need for several specific fea- 
tures previously considered in the literature for related 
networks. As already mentioned, oscillations require nei- 
ther an additional positive feedback loop |22|. nor self- 
activation of gene A |23j nor a highly cooperative binding 
|24j of transcription factors to gene promoters. Similarly, 
an explicit delay is not needed although delays of various 
origins have been considered in models of oscillatory ge- 
netic networks in different contexts, as further discussed 
below. 

Besides the motif topology, one feature of the present 
MFL model that renders this possible is the explicit de- 
scription of transcription and translation. The condensed 
representation of protein production by a single effective 
process that is often used, would necessitate supplemen- 
tary interactions to make oscillations possible. The de- 
scription of mRNA was for instance omitted in the first 
version of a previously proposed evolutionary algorithm 
5j. Oscillatory networks with a core MFL motifs were 
nonetheless produced but much less easily than in a re- 
fined second version that incorporates mRNA dynamics. 
This role of mRNA dynamics may well be worth bear- 
ing in mind when considering the elaborate regulation of 
translation that is presently being uncovered |25j . 
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B. Biological MFL oscillators and switches. 

Given the proposed functions for the MFL motif, it is 
of primary interest to know whether cases of its use as an 
oscillator or a switch are already documented. The os- 
cillator function seems the clearest. Circadian clocks are 
genetic oscillators which generate endogenous rhythms 
with a period close to 24 hours and which are locked to an 
exact 24 hours period by the alternation of day and night. 
A MFL motif is found at the core of all known eucaryotic 
circadian oscillators. Taking for illustrative purposes the 
fungus Neurospora crassa, one of the best studied model 
organisms, it has been established that a dimer of White- 
Collar proteins (WCC) plays the role of A and activates 
the transcription of the frequency gene, the gene in 
this example. In turn, the FRQ protein interacts with 
WCC and prevents this activation [26| . Analogous motifs 
are found in the circadian genetic networks of flies and 
mammals, as shown in Table [I] The p53/Mdm2 mod- 
ule [23 provides a different example. The p53 protein 
is a key tumor suppressor protein and an important role 
is played in its regulation by the Mdm2 protein. On the 
one hand, p53, similarly to the A protein, binds to Mdm2 
gene an activates its transcription. On the other hand, 
the Mdm2 protein, as the B protein, binds to p53 and 
both blocks transcriptional activation by p53 and pro- 
motes its rapid proteolytic degradation. Cells exposed 
to stress have indeed been observed to present oscilla- 
tions in both p53 and Mdm2 levels 27]. In both these 
examples, many other genetic interactions exist besides 
the above described MFL motifs and have been thought 
to provide delays necessary for oscillations. These have 
been postulated to come from chains of phosphorylations 
in circadian networks |28j or from a yet unknown inter- 
mediate species [2^ in the p53-Mdm2 system. The re- 
alization that the MFL can lead to oscillations without 
further interactions may be useful in suggesting alter- 
native models of these particular genetic networks or in 
reassessing the role of known interactions. For instance, 
this has led one of us to formulate a model of the Neu- 
rospora crassa circadian clock |29l ] in which kinases and 
phosphorylation influence the cycle period by modifying 
proteins degradation rates but are not needed to create 
key delays. The prevalence of the MFL motif probably 
arises from the usefulness of negative feedback but it may 
also imply that oscillations in genetic network are more 
common than usually thought. 

The evidence for uses of the MFL motif as a switch 
appears less clear-cut at present. The classic case of the 
E. Coli lactose operon [3(j can be thought as an effective 
version of a bistable MFL. The lac repressor represses the 
lac gene and plays the role of A in the MFL. The lac gene 
directs the production of a membrane permease which 
itself drives the absorption of external lactose. Since al- 
lolactose binding to the lac repressor blocks its transcrip- 
tional activity, allolactose effectively plays the same role 
as B in the MFL. The search for a more direct switch 
example has led us to consider the different MFL motifs 



A Species 


B Species 


Main biological functions 


Bistable systems 


Lac Repressor 


Allolactose 


Lactose metabolism 


Oscillators 


WC-1, WC-2 


FRQ 


Neurospora circadian clock 


dCLK 


PER, TIM 


Drospophila circadian clock 


CLOCK, BMAL 


PER, CRY 


Mammals circadian clock 


p53 


Mdm2 


Stress response oscillators 


Detected Yeast MFLs 


Swi6 


Swi4 


Gl/S transition 


Gal4 


Gal80 


Galactose metabolism 


Gal4 


Gal3 


Galactose metabolism 


Gal4 


Gall 


Galactose metabolism 


Imel 


Rim 11 


Meiose activation 


Ume6 


Imel 


Meiose activation 


Stel2 


Fus3 


Pseudohyphal growth 


Stel2 


Farl 


Pseudohyphal growth 


Cbfl 


Met28 


Sulfur metabolism 


Met4 


Met28 


Sulfur metabolism 


Swi4 


Clb2 


Cell cycle 


Mbpl 


Clb5 


Gl/S transition 


Stbl 


Clnl 




Stbl 


Cln2 





TABLE I: Some biological examples of MFL motifs. The yeast 
motifs are reproduced from ref. Q ■ An annotation taken from 
the SGD database [33] has been added. 



in yeast as reported in Q and reproduced in Table [I] It 
was noted in Q that most of these MFL motifs were cen- 
tral modules in biochemical pathways. Most of them take 
part in the cell cycle or in differentiation pathways, which 
certainly are processes where the cell switches from one 
function to another. However, we have found it difficult 
to disentangle the MFL's from numerous other known in- 
teractions and to confidently conclude that any of these 
motifs implements the proposed switch function. Simi- 
lar difficulties have been recently emphasized for general 
motifs determined on purely statistical grounds and have 
led to question their functional significance |3l|. In the 
present case, the difficulty of identifying biological cases 
of MFL switches, of course arises from our own very par- 
tial knowledge of the networks listed in Table I, but it 
may also be due to the fact that the possible role of the 
MFL module as a switch was not fully realized in previ- 
ous investigations. The present study will hopefully help 
and trigger direct experimental investigations of these 
questions. 



APPENDIX A: TRANSITIONS BETWEEN THE 
TWO PHASES OF AN OSCILLATION CYCLE 
AND VS CORRECTION TO THE PERIOD. 

In the following, we analyze the transitions between 
phase I and II of an oscillation period and use matched 
asymptotics to precisely justify the assumptions made 
in section Ivl This also provides the leading order (V~8) 
correction to the zeroth order results of section Ivl 
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1. From phase I to phase II. 

We consider the transition between the end of the high 
A/low B phase I and the start of the high B/low A phase 
II. This protein switch occurs on a time scale of order 
S/S r . On this fast time scale the mRNA concentration 
does not have time to change. Thus, in Eqs. 1( . we 
assume r — r 2 and introduce r — t/S. Eqs. I(8I9|) simply 
become at dominant order 



= r 2 -AB 
= 1-AB 
The imposed boundary conditions are 



dB_ 

dT 
dA 

"d7 



(Al) 
(A2) 



A (l-r 2 )(r-n)+o(l), B 
B -> (r 2 -l)(T-T 2 )+o(l), A ■ 



at r = -oo (A3) 
at t — > +oo (A4) 



with the dimensionless transcript concentration r 2 > 1 
(see section IV B 1(1 and n and t 2 two constants to be 
determined. Eqs. (|A1IA2|I are autonomous in time and 
therefore invariant by time translation. For general non 
integrable equations, numerical integration would be re- 
quired to obtain the difference r 2 — t%. Here, however, 
the difference of protein concentrations, A — B, is easily 
integrated across the transition region. Fixing the time 
origin at the instant when A — B, one readily obtains 
the exact formula 



A-B=(l-r 2 ) 



(A5) 



Comparison with Eqs. (|A3IA4() shows that with this 
choice of time origin t% = r 2 = (and more generally 
Ti = r 2 ). The asymptotic conditions l|A3IA4|) coincide 
with the limiting behavior of a / (t) / 8 at the end of phase 
I near t — t\, and with the limiting behavior of bn(t) 
at the beginning of phase II near t = 0, provided that 
o>i{t\ ) = bn(p) = as was required in the main text. 
Thus, Eqs. (|A1IIA2|) provide a uniform approximation of 
A(t) and B(t) throughout the transition from phase I to 
phase II. Given the exact result (|A5|> . the integration of 
Eqs. (| All IA2|) can in fact be replaced to the integration 
of the single following Riccati equation: 



^ = 1 - [A + (ra - l)r]A 
dr 



(A6) 



A comparison between Eq. (|A6|) and full numerical evo- 
lution is shown in Fig. 1101 There is one subtlety, however. 
The quasi-equilibrium approximation for g [Eq. (|42|l ] di- 
verges at the end of phase I when A — > and g becomes 
larger than 1 before the transition region (|A1IA2|I of or- 
der 5. This clearly signals the breakdown of the approxi- 
mation (|42|l in a larger intermediate region at the end of 
phase I. In a region of size t\ — t <~ Vo, the evolution of 
g reduces to 



M =e ( 1 - g 5A Q 



(A7) 



Indeed, taking a = (1 — r 2 )(t\ —t) + o(V~5) shows that all 
three terms kept in Eq. QA7JI are of the same magnitude 
when g ~ t\ — t ~ VS and that the additional term 
8g in Eq. © is negligible. Explicit integration in this 
intermediate region gives 

g(ti) — g(t) exp [— n 2 (t\ — f) 2 ] = 6 duexp(n 2 u 2 ) 



(A8) 

where k 2 = 6(1 — r 2 )/(2AoS). In particular, the limit 
t — * — oo, determines the value of g at the end of phase I 
and beginning of phase II 



itA q 65 



92 



2(r 2 - 1) 



(A9) 



Eq. I|A9|) is the source of a correction of order y/~5 to 
the asymptotic result l|55|) as will be shown below. 



2. From phase II to phase I. 

The transition from phase II to phase I can be ana- 
lyzed quite similarly to the transition from phase I to 
phase II. In a time of order S/S r at the end of phase II, 
r has no time to change and the transition is described 
by Eqs. I|A1IA2(I with r 2 replaced by < 1 (see section 
IV B 1J) and the boundary conditions 



B 
A 



(n-i)(r- 

(l-ri)(r- 



t[) + o{1), A 
7^) + o(l), B 



at r : 
at r 



-oo (A10) 
-oo (All) 



Again, r[ 



if the time origin is taken when 



A = B. This transition regime is followed by a longer 
transition on a time scale t — £ 2 ~ >/5 where g decreases 
from a value of order one to values of order S character- 
istic of phase I. This reflects the fact that the binding of 
protein A to the promoter of gene gb is not instantaneous 
and requires some accumulation of protein A at the be- 
ginning of phase I. Taking A = (1 — ri)t/5, g evolution 
is described in this intermediate regime by 



dg 
dt 



l-t(l 



The corresponding evolution of g is 



g(t) = exp (—Kit 2 ) 



gi + 6 du exp (ki« 2 ) 
Jo 



(A12) 



(A13) 



where k± = 6(1 — ri)/(2AoS). It indeed describes the 
transition from g\ = 1 — exp(#£ 2 ), the value of g at the 
end of phase II, to the quasi-equilibrium regime since the 
asymptotic behavior of Eq. (|A13|I for t ^J~8 is 



9(t) 



And 



(i-n)t 



(A14) 
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3. Dominant correction to the oscillation period. 

The form of Eqs. (l(iF)ll could lead to think that the 
oscillation period has an expansion in powers of the small 
parameter S. However, here as often, boundary layers 
lead to more complicated expansions. The two transition 
regimes of duration y/5 give rise to a dominant correction 
of order y/~8 to the oscillation period. 

As pointed out above, due to the transition regime at 
the end of phase I, the starting value of g in phase II is 
g 2 [Eq. l|A9|l ] of order y/S. Therefore, the first Eq. 
is replaced at order y/~5 by the corrected evolution, 

9ii At) = 1 - (1 - 92) exp(-9t) (A15) 

This leads in turn to corrected evolutions ru iC (t) an d 
bii. c (t) for the mRNA and B protein concentrations that 
are obtained by simply replacing (pi — po) by (p% — po)(l — 
g 2 ) in the expressions l|52[l of the zeroth-order expressions 
r H (t) and b H (t). 

The other transition regime between phase II and 
phase I does not lead to modifications of the form of 
Eq. I|46I47J) describing the evolutions of mRNA and A 
protein during the bulk of phase I. Instead, it results in 
a difference of order y/S between rn^fa) and r\ the ef- 
fective concentrations of transcripts at the start of phase 
I. The transcript concentration depends on the promoter 
states at previous times. So, integrating Eq. (JJJ) from the 
end of phase II, one obtains 

r(t) = pi + [rnAte) - Pi] exp(-t) 

+ (Po - Pi)cxp(-i) duexp(u)g(u) (A16) 
Jo 



Replacing g(t) by its expression (|A12(1 during the tran- 
sition regime at the start of phase I leads for t y/8 to 



r{t) = Pi + [rn, c (h) - 91 J -. — {pi - pa) - Pi] exp(-t) 

V 4k 1 

(A17) 

Namely, the previous Eq. l(4^|) but with 



ri=rn iC (t 2 )- giJ — (pi- po) (A18) 

Note that only the gi-term between the square brackets 
in Eq. (|A12|I contribute to Eq. (IA17f) . the integral term 
giving a subdominant contribution. 

Finally, the durations T\ and T 2 of both phases and the 
total period T r = T\ + T 2 are obtained to yfb accuracy 
by solving as before 

ri(t 1 )=r 2 , a/(ti) = (A19) 

together with the corrected version of Eq.JMJ), namely 
Eq. PT%|) and b TI , c (t 2 ) = . 

The asymptotic description of gb evolution agrees well 
with numerical simulations of the full MFL model as 
shown in Fig. 1111 The y/S correction terms lengthen the 
zeroth-order asymptotic estimation. For moderate value 
of (5, a numerically comparable contribution is however 
coming from higher-order terms that tend to lengthen 
the period [coming, for instance, from the breakdown of 
the quasi-equilibrium approximation for gb Eq. (|42H ] as 
can be seen in Fig. |S] 
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FIG. 2: The proposed model of the MFL motif. The Greek 
letters denote the model different kinetic constants. The 
large crosses symbolize the degradation of the corresponding 
species. 



FIG. 1: Schematic representation of the over-represented 
MFL motif |4|. The bold arrow represents a transcriptional 
regulation interaction and the dashed double arrow a protein- 
protein interaction. 
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FIG. 3: Different dynamical regimes of the MFL as po and 
pi are varied. The borders between different regimes are 
shown as computed from numerical solutions of Eq. 11811 and 
(ED for 6» = 1.33 (full lines) or 9 = 26.6 (d ashed) as well 
as given by the asymptotics expansions [Eqs. I28I30I36I 1411 1 
for 9 — 1.33 (+ symbols). Here, and in all following fig- 
ures, except when explicitly specified, the others parameters 
are [A]o = 4Qmol,f3 = 3 min" 1 , 6 r — 0.03 miri~ , 8a = 
0.01 min -1 , 6b = 0.01 min" 1 , pA = 100 rn.ol.min~ 1 , 
7 = 1 mol~ x .rain~ x where mol stands for molecules and 
mill for minutes. The corresponding dimensionless parame- 
ters are 6 = 0.003, d a = d b = 0.33, A = 4, p = 0.31. 
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FIG. 4: In the bistable regime a graded increase of pa, the 
production of A, results in a jump in A concentration. The 
parameters are p/ — 0.2 mol.min~ pt = 2 mol.mivT 1 , 
and 6 = 0.04 min~ . The other parameters are as in Fig.[!|] 
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FIG. 5: The oscillation domain for the same parameters as in 
Fig.[3]but f or a larger domain of pi, the ratio of activated pro- 
duction of B to that of A. As in Fig. the short-dashed line 
marks the boundary of the domain where the steady stable is 
unstable to oscillations for 9 = 26.6. Numerical simulations 
(diamonds) show that the oscillating regime is stable up to 
the long-dashed line. Therefore, both the limit-cycle and the 
fixed point are stable attractors for 9 = 26.6 in the region 
between the short-dashed and long-dashed lines. 
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FIG. 6: MFL in the oscillatory regime. The concentrations of 
the different species are shown as a function of time and dis- 
play sustained oscillations. Constants are the same as in Fig. 
^except that pa = 100 mol.min' 1 , pj = 0.1 mol.min~ 
pb = 5 mol.min^ 1 . 
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FIG. 7: Distinction between the two phases for the adimen- 
sioned equations [Eqs. (161911 1 . Top panel : oscillation of r. 
Bottom panel : oscillations for A (full line) and B (dotted). 
In order to clearly depict phase I and II of an oscillation cycle, 
the parameters are here chosen so that the two phases have 
similar durations : 8 = 0.001, pi = 1.5,po = 0, 6 = 2, Aq = 1, 
fi = 1, d a = db = 0. 




FIG. 8: Comparison between the lowest order theoretical 
rescaled period T r (bold line), the approximate expression for 
large pi [Eq. 1611 1 (dashed line) and numerically computed 
periods for different values of 8, as a function of parameter pi . 
Other parameters are 9 — 10, Aq = 10, jit = 1, d a = db = 0. 
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FIG. 9: Comparison between the full dynamics (+) sym- 
bols and the asymptotic description at order yS (dashed and 
dotted lines) during one oscillation cycle, (a) Concentration 
of mRNA behavior . Phase I asymptotics (dotted lines) and 
phase II asymptotics (dashed lines) are shown separately. The 
product A.B//3 is also plotted (full line). The quasi-static as- 
sumption is seen to be well-satisfied away from the transition 
regions between the two phases, (b) Concentration of protein 
B. (c) Concentration of protein A. Parameters are as in Fig. 
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FIG. 10: The evolution of A given by the Riccati Eq. ijUjl 
(+) is compared to that given by the complete MFL dynam- 
ics (dashed) during the transition from phase I to phase II. 
Parameters are as in Fig. [(J 
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FIG. 11: Comparison of [gt] obtained from a full numerical 
integration (+) with the uniform approximation obtained by 
matching the different transition regimes (bold line). The 
exponential relaxation in phase II [Eq. (IA15H is also shown 
(dashed). An enlargment of phase I and the transition regimes 
is shown in the inset. Parameters are as in Fig. |S| 



